function eq = solve_eq_ss(param,glob,options)

% Finds initial steady state
% With consumption normalized to 1

 D1 = 0;
 D2 = 10;
 Dmax = D2;
 
% Impose that C = 1. 
param.C = 1;   

Ds    = nan(options.itermax_DC,1);

if nargin == 5
    param.D = eq0.D_out;
    options.cresult = eq0.c;
else
    param.D = .459;
end
    
param.w = 1;

tic   

% Use bisection to find the value of D such that vE = cE
tosolve = @(D) solve_cL(D, param,glob,options);
param.D = bisect(tosolve, 1e-3, 3);
[~,eq] = solve_cL(param.D, param,glob,options);

% Now find the mass of entrants that is consistent with that value of D
eq = setup_transition_matrices(eq, param,glob,options); % set up transition matrices
tosolve = @(M) find_stationary_dist(M, eq, param,glob,options);
M = bisect(tosolve, 1e-3, 100);

[diff,eq] = find_stationary_dist(M, eq, param,glob,options);
toc

    subplot(3,1,3)
    plot(glob.bgridf, eq.Lp)
    title('Distribution of labor'); drawnow;
    options.cresult = eq.c;
 
   
%    param.D = options.damp_DC*param.D + (1-options.damp_DC)*eq.D_out;
    eq.D = param.D;
    eq.M = M;
end